library(TTR)
library(quantmod)
library(ggplot2); library(readr); library(dplyr); library(tidyverse);
library(gridExtra);library(ggExtra); library(knitr);
library(caret); library(gridExtra); library(lubridate); library(stringr);
library(PerformanceAnalytics); library(xts); library(tidyquant); library(tidyverse); library(data.table); library(chron)
us_mass <- read.csv("US_Mass_Shootings_May_24_2022.csv")
#Select those variables that are needed for our analysis
us_mass02 <- us_mass %>% select(location, date, fatalities, injured, total_victims)
str(us_mass02)
## 'data.frame': 128 obs. of 5 variables:
## $ location : chr "Uvalde, Texas" "Buffalo, New York" "Sacramento, California" "Oxford, Michigan" ...
## $ date : chr "5/24/22" "5/14/22" "2/28/22" "11/30/21" ...
## $ fatalities : int 15 10 4 4 9 8 4 10 8 4 ...
## $ injured : chr "-" "3" "0" "7" ...
## $ total_victims: chr "-" "13" "4" "11" ...
# change date format to standardize dates
us_mass02 <- us_mass02 %>% mutate(date = mdy(date)) %>%
#Shooting was/is in a holiday date is move to next date
mutate(date = ifelse(is.holiday(date), date+1,date)) %>%
# Moving any dates to the following Monday
mutate(weekday=wday(date)) %>%
mutate(date1 = ifelse(weekday==1,(date+1), ifelse(weekday==7, (date+2), (date)))) %>%
mutate(date = ifelse(is.holiday(date), date+1,date)) %>%
mutate(date = as.Date(date1)) %>%
mutate(year = year(date)) %>%
mutate(state = ifelse(location=="Washington, D.C.", "Washington, D.C.",str_extract(location,
"[\\s]*[A-Z][a-z]+[\\s]*[[A-Z][a-z]+]*$"))) %>%
mutate(fatalities = as.integer(fatalities)) %>%
mutate(injured = as.integer(injured)) %>%
mutate(total_victims = as.integer(total_victims)) %>%
select(location, state, date, year, fatalities, injured, total_victims)
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
## Update Uvalde, Texas data
us_mass02$fatalities[us_mass02$location=="Uvalde, Texas"] <- as.integer(22)
us_mass02$injured[us_mass02$location=="Uvalde, Texas"] <- as.integer(17)
us_mass02$total_victims[us_mass02$location=="Uvalde, Texas"] <- as.integer(39)
## Add Highland Park, Illinois mass shooting
us_mass02 <- us_mass02 %>% add_row(location = "Highland Park, Illinois", state = "Illinois", date = mdy(07042022), year= as.integer(2022), fatalities = as.integer(8), injured=as.integer(29), total_victims=as.integer(37), .before=1)
summary(us_mass02)
## location state date year
## Length:129 Length:129 Min. :1982-08-20 Min. :1982
## Class :character Class :character 1st Qu.:2001-02-05 1st Qu.:2001
## Mode :character Mode :character Median :2013-06-07 Median :2013
## Mean :2009-10-08 Mean :2009
## 3rd Qu.:2018-02-14 3rd Qu.:2018
## Max. :2022-07-04 Max. :2022
## fatalities injured total_victims
## Min. : 3.000 Min. : 0.00 Min. : 3.00
## 1st Qu.: 4.000 1st Qu.: 1.00 1st Qu.: 7.00
## Median : 6.000 Median : 3.00 Median : 10.00
## Mean : 8.093 Mean : 11.69 Mean : 19.78
## 3rd Qu.: 9.000 3rd Qu.: 10.00 3rd Qu.: 18.00
## Max. :58.000 Max. :546.00 Max. :604.00
Group mass shooting data by state and year
usmass_yearly <- us_mass02 %>% group_by(year) %>%
summarise( fatalities=sum(fatalities, na.rm = TRUE), injured=sum(injured, na.rm = TRUE), total_victims=sum(total_victims, na.rm = TRUE))
shootings_yrl <- us_mass02 %>% count(year)
colnames(shootings_yrl) <- c("year", "mass_shootings")
usmass_yearly <- merge(usmass_yearly, shootings_yrl, by="year")
usmass_yearly$intensity <- usmass_yearly$total_victims/usmass_yearly$mass_shootings
tail(usmass_yearly)
## year fatalities injured total_victims mass_shootings intensity
## 33 2017 117 587 704 11 64.000000
## 34 2018 80 70 150 12 12.500000
## 35 2019 73 112 185 10 18.500000
## 36 2020 9 0 9 2 4.500000
## 37 2021 43 16 59 6 9.833333
## 38 2022 44 49 93 4 23.250000
usmass_state <- us_mass02 %>% group_by(state) %>%
summarise( fatalities=sum(fatalities, na.rm = TRUE), injured=sum(injured, na.rm = TRUE), total_victims=sum(total_victims, na.rm = TRUE))
shootings_state <- us_mass02 %>% count(state)
colnames(shootings_state) <- c("state", "mass_shootings")
usmass_state <- merge(usmass_state, shootings_state, by="state")
usmass_state <- usmass_state %>% arrange(desc(mass_shootings))
head(usmass_state)
## state fatalities injured total_victims mass_shootings
## 1 California 157 160 317 23
## 2 Florida 126 109 235 12
## 3 Texas 152 183 335 12
## 4 Colorado 48 104 152 7
## 5 Washington 37 28 65 7
## 6 New York 40 28 68 5
Plot of total fatalities grouped by state
subset_usmass_state1 <- usmass_state %>% arrange(desc(total_victims)) %>% slice(1:12)
temp_masssate <- subset_usmass_state1[,c("state", "injured", "fatalities")] %>% reshape2::melt()
## Using state as id variables
usshoot_state_plot <- ggplot(data = temp_masssate, aes(y=reorder(state, (value)), x=value, fill=variable))
usshoot_state_plot <- usshoot_state_plot + geom_col(position="stack", stat="identity")
## Warning: Ignoring unknown parameters: stat
usshoot_state_plot <- usshoot_state_plot + scale_x_continuous(breaks = seq(0, 700, by = 100))
usshoot_state_plot <- usshoot_state_plot + labs(title="US Mass Shooting Total Victims per State",
subtitle = "Kaggle US Mass Shootings",
caption = "Team 89")
usshoot_state_plot <- usshoot_state_plot + ylab("State") + xlab("Total Victims")
usshoot_state_plot <- usshoot_state_plot + theme(plot.title = element_text(size = 15))
usshoot_state_plot <- usshoot_state_plot + theme(axis.title = element_text(size = 15))
usshoot_state_plot <- usshoot_state_plot + theme(axis.text = element_text(size = 15))
usshoot_state_plot <- usshoot_state_plot + theme(legend.title = element_blank())
usshoot_state_plot <- usshoot_state_plot + theme(legend.position = "right")
usshoot_state_plot
ggsave("89plot01.png")
## Saving 7 x 5 in image
Plot of number of mass shootings grouped by state
#temp_masssate <- usmass_state[,c("state", "injured", "fatalities")] %>% reshape2::melt()
subset_usmass_state2 <- usmass_state %>% arrange(desc(mass_shootings)) %>% slice(1:12)
usshoot_state_plot1 <- ggplot(data = subset_usmass_state2, aes(y=reorder(state, mass_shootings), x=mass_shootings))
usshoot_state_plot1 <- usshoot_state_plot1 + geom_col(position="stack", stat="identity", fill="red")
## Warning: Ignoring unknown parameters: stat
usshoot_state_plot1 <- usshoot_state_plot1 + scale_x_continuous(breaks = seq(0, 30, by = 5))
usshoot_state_plot1 <- usshoot_state_plot1 + labs(title="US Total of Mass Shootings per State",
subtitle = "Kaggle US Mass Shootings",
caption = "Team 89")
usshoot_state_plot1 <- usshoot_state_plot1 + ylab("State") + xlab("Total Mass Shootings")
usshoot_state_plot1 <- usshoot_state_plot1 + theme(plot.title = element_text(size = 15))
usshoot_state_plot1 <- usshoot_state_plot1 + theme(axis.title = element_text(size = 15))
usshoot_state_plot1 <- usshoot_state_plot1 + theme(axis.text = element_text(size = 15))
usshoot_state_plot1 <- usshoot_state_plot1 + theme(legend.title = element_blank())
usshoot_state_plot1 <- usshoot_state_plot1 + theme(legend.position = "right")
usshoot_state_plot1
ggsave("89plot02.png")
## Saving 7 x 5 in image
Plot of number of mass shooting per year
usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=mass_shootings), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=mass_shootings), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Annual Mass Shootings",
subtitle = "Kaggle US Mass Shootings",
caption = "Team 89")
usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'
ggsave("89plot03.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
Plot of annual fatalities due to mass shootings
usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=fatalities), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=fatalities), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Annual Mass Shooting Fatalities",
subtitle = "Kaggle US Mass Shootings",
caption = "Team 89")
usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'
ggsave("89plot04.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
Number of victims per shooting on function of year
usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=intensity), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=intensity), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Number of Victims per mass shooting",
subtitle = "Kaggle US Mass Shootings",
caption = "Team 89")
usshoot_plot1 <- usshoot_plot1 + ylab("Victims per mass shooting")
usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'
ggsave("89plot05.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
Obtain stock history for the companies that are been evaluated using tq_get function
date_str <- "2009-10-01"
#stocks <- c("swbi")
stocks <- c("rgr", "swbi")
#stocks <- c("rgr", "swbi", "oln", "npk", "bgfv", "axon" )
gun_co_stocks <- stocks %>% tq_get(get="stock.prices", from=date_str, to=Sys.Date())
head(gun_co_stocks)
## # A tibble: 6 × 8
## symbol date open high low close volume adjusted
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 rgr 2009-10-01 12.9 13.0 12.4 12.5 252500 8.38
## 2 rgr 2009-10-02 12.4 12.7 11.8 12.3 254000 8.25
## 3 rgr 2009-10-05 12.3 12.5 12.2 12.2 120300 8.20
## 4 rgr 2009-10-06 12.2 12.8 12.2 12.5 111300 8.41
## 5 rgr 2009-10-07 12.5 12.8 12.4 12.6 99300 8.45
## 6 rgr 2009-10-08 12.7 13.0 12.5 12.8 144400 8.60
The historical stock price for analyzed public companies, since October 2009
gun_stocks_plot <- ggplot(data = gun_co_stocks, aes(x=date, y=adjusted, group=symbol))
gun_stocks_plot <- gun_stocks_plot + geom_line(aes(color=symbol))
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom")
gun_stocks_plot <- gun_stocks_plot + labs(title = "Adjusted Stock Price US Gun Manufactures",
subtitle = "Source: Yahoo Finance",
caption = "team 89")
gun_stocks_plot <- gun_stocks_plot + ylab("Adjusted Stock Price")
gun_stocks_plot <- gun_stocks_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_stocks_plot
ggsave("89plot06B.png")
## Saving 7 x 5 in image
Generate a vector with dates of mass shootings so we can use it in the following exploratory plots
us_mass03 <- us_mass02 %>% filter(date>="2010-01-01" & total_victims>=7)
us_mass_dates <- us_mass03[,3]
Generate plots including adjusted stock price and dates of mass shootings (vertical dash lines)
gun_stocks_plot <- ggplot(data = gun_co_stocks, aes(x=date, y=adjusted, group=symbol))
gun_stocks_plot <- gun_stocks_plot + geom_line(aes(color=symbol), size=0.75)
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom")
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom", legend.key.width=unit(0.1,"npc"))
## Vertical lines to show mass shootings
gun_stocks_plot <- gun_stocks_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_stocks_plot <- gun_stocks_plot + labs(title = "Adjusted Stock Price US Gun Manufactures",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_stocks_plot <- gun_stocks_plot + ylab("Adjusted Stock Price")
gun_stocks_plot <- gun_stocks_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_stocks_plot
ggsave("89plot07B.png")
## Saving 7 x 5 in image
Obtain the returns for each of companies using tq_transmute function from tidyquant package
gun_co_returns <- gun_co_stocks %>% group_by(symbol) %>% tq_transmute(select = adjusted,
mutate_fun = periodReturn,
period="daily",
col_rename = "Ra")
head(gun_co_returns)
## # A tibble: 6 × 3
## # Groups: symbol [1]
## symbol date Ra
## <chr> <date> <dbl>
## 1 rgr 2009-10-01 0
## 2 rgr 2009-10-02 -0.0160
## 3 rgr 2009-10-05 -0.00570
## 4 rgr 2009-10-06 0.0254
## 5 rgr 2009-10-07 0.00559
## 6 rgr 2009-10-08 0.0175
Generate a plot showing return of each of the companies and dates of mass shootings (vertical dash lines)
gun_returns_plot <- ggplot(data = gun_co_returns, aes(x=date, y=Ra, group=symbol))
gun_returns_plot <- gun_returns_plot + geom_line(aes(color=symbol))
gun_returns_plot <- gun_returns_plot + theme(legend.position = "bottom")
gun_returns_plot <- gun_returns_plot + theme(legend.position = "bottom", legend.key.width=unit(0.1,"npc"))
## Vertical lines to show mass shootings
gun_returns_plot <- gun_returns_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_returns_plot <- gun_returns_plot + labs(title = "Returns of US Gun Manufactures",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_returns_plot <- gun_returns_plot + ylab("Ra") + ylim(-0.2, 0.2)
gun_returns_plot <- gun_returns_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_returns_plot
ggsave("89plot08B.png")
## Saving 7 x 5 in image
rvest https://cran.r-project.org/web/packages/rvest/rvest.pdf package is used to harvest (scrap) summary page of each company in Yahoo Finance. Once the market values are obtain a weights are determine to set up portafolio
library(rvest)
## Warning: package 'rvest' was built under R version 4.2.1
##
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
##
## guess_encoding
# Define stock name
#stocks <- c("rgr", "swbi")
#stocks <- c("rgr", "swbi", "oln", "npk", "bgfv", "axon" )
Mkt_Caps <- c()
for(s in stocks){
# Extract and transform data
temp_Co_Metrics <- paste0("https://finance.yahoo.com/quote/", s) %>%
read_html() %>% html_table() %>% map_df(bind_cols) %>%
t() %>% as_tibble()
temp_Mkt_Cap <- temp_Co_Metrics[2,9]
temp_Cap_Unit <- str_extract(temp_Mkt_Cap, "[A-Z]")
temp_Cap_Amount <- as.numeric(str_extract(temp_Mkt_Cap, "\\d*\\.\\d*"))
temp_Cap_Value <- case_when(temp_Cap_Unit== "T" ~ temp_Cap_Amount*1000*1000,
temp_Cap_Unit== "B" ~ temp_Cap_Amount*1000,
temp_Cap_Unit=="M" ~ temp_Cap_Amount*1,
is.na(temp_Cap_Unit) ~ temp_Cap_Amount/(1000*1000))
Mkt_Caps <- append(Mkt_Caps, temp_Cap_Value)
}
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
gun_wts <- Mkt_Caps/sum(Mkt_Caps)
#gun_wts <- c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)
co_summary <- data.frame(stocks, Mkt_Caps, gun_wts)
co_summary
## stocks Mkt_Caps gun_wts
## 1 rgr 1112.000 0.6413195
## 2 swbi 621.925 0.3586805
A data set for the portafolio was created as follows:
gun_co_portfolio_returns <- gun_co_returns %>%
tq_portfolio(assets_col = symbol,
returns_col = Ra,
weights = gun_wts,
col_rename = "Ra")
## Warning: `spread_()` was deprecated in tidyr 1.2.0.
## Please use `spread()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% tq_mutate(select = Ra, mutate_fun = SMA, n=2) %>%
rename(SMA.2day = SMA) %>%
tq_mutate(select = Ra, mutate_fun = SMA, n=15) %>%
rename(SMA.15day = SMA)
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(gunDD = Drawdowns(Ra)) %>%
mutate(gunDD=as.numeric(gunDD)) %>%
mutate(gunDD_2day=lead(gunDD,2)) %>%
mutate(gunDD_3day=lead(gunDD,3)) %>%
mutate(gunDD_5day=lead(gunDD,5)) %>%
mutate(gunDD_10day=lead(gunDD,10)) %>%
mutate(gunDD_15day=lead(gunDD,15)) %>%
mutate(gunDD_20day=lead(gunDD,20))
## Warning in merge.zoo(fx, .xts(, .index(x))): Index vectors are of different
## classes: integer POSIXct
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(Ra_2day=frollmean(Ra,2)) %>%
mutate(Ra_3day=frollmean(Ra,3)) %>%
mutate(Ra_5day=frollmean(Ra,5)) %>%
mutate(Ra_10day=frollmean(Ra,10))
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(Mom_2day=lead(Ra,2) - Ra) %>%
mutate(Mom_3day=lead(Ra,3) - Ra) %>%
mutate(Mom_5day=lead(Ra,5) - Ra) %>%
mutate(Mom_10day=lead(Ra,10) - Ra) %>%
mutate(Mom_15day=lead(Ra,15) - Ra) %>%
mutate(Mom_20day=lead(Ra,20) - Ra)
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% filter(date>="2010-01-01")
head(gun_co_portfolio_returns,20)
## # A tibble: 20 × 21
## date Ra SMA.2day SMA.15day gunDD gunDD_2day gunDD_3day
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-04 0.0398 0.00672 0.00234 -0.209 -0.174 -0.173
## 2 2010-01-05 0.0255 0.0326 0.00390 -0.189 -0.173 -0.179
## 3 2010-01-06 0.0177 0.0216 0.00502 -0.174 -0.179 -0.177
## 4 2010-01-07 0.00170 0.00969 0.00560 -0.173 -0.177 -0.187
## 5 2010-01-08 -0.00690 -0.00260 0.00413 -0.179 -0.187 -0.182
## 6 2010-01-11 0.00154 -0.00268 0.00489 -0.177 -0.182 -0.181
## 7 2010-01-12 -0.0123 -0.00539 0.00340 -0.187 -0.181 -0.157
## 8 2010-01-13 0.00697 -0.00268 0.00346 -0.182 -0.157 -0.141
## 9 2010-01-14 0.000996 0.00398 0.00385 -0.181 -0.141 -0.164
## 10 2010-01-15 0.0294 0.0152 0.00590 -0.157 -0.164 -0.178
## 11 2010-01-19 0.0186 0.0240 0.00720 -0.141 -0.178 -0.183
## 12 2010-01-20 -0.0270 -0.00421 0.00633 -0.164 -0.183 -0.180
## 13 2010-01-21 -0.0159 -0.0215 0.00484 -0.178 -0.180 -0.181
## 14 2010-01-22 -0.00617 -0.0110 0.00317 -0.183 -0.181 -0.181
## 15 2010-01-25 0.00366 -0.00125 0.00517 -0.180 -0.181 -0.193
## 16 2010-01-26 -0.00193 0.000863 0.00239 -0.181 -0.193 -0.212
## 17 2010-01-27 0.000884 -0.000523 0.000749 -0.181 -0.212 -0.197
## 18 2010-01-28 -0.0145 -0.00683 -0.00140 -0.193 -0.197 -0.192
## 19 2010-01-29 -0.0242 -0.0194 -0.00312 -0.212 -0.192 -0.184
## 20 2010-02-01 0.0194 -0.00237 -0.00137 -0.197 -0.184 -0.207
## # … with 14 more variables: gunDD_5day <dbl>, gunDD_10day <dbl>,
## # gunDD_15day <dbl>, gunDD_20day <dbl>, Ra_2day <dbl>, Ra_3day <dbl>,
## # Ra_5day <dbl>, Ra_10day <dbl>, Mom_2day <dbl>, Mom_3day <dbl>,
## # Mom_5day <dbl>, Mom_10day <dbl>, Mom_15day <dbl>, Mom_20day <dbl>
Gun portfolio statistics using table.Stats from PerformaceAnalytics package https://www.rdocumentation.org/packages/PerformanceAnalytics/versions/2.0.4/topics/table.Stats
gun_portfolio_stats <- table.Stats(gun_co_portfolio_returns$Ra)
gun_portfolio_stats
##
## Observations 3160.0000
## NAs 0.0000
## Minimum -0.1771
## Quartile 1 -0.0107
## Median 0.0012
## Arithmetic Mean 0.0009
## Geometric Mean 0.0006
## Quartile 3 0.0125
## Maximum 0.1397
## SE Mean 0.0004
## LCL Mean (0.95) 0.0001
## UCL Mean (0.95) 0.0017
## Variance 0.0006
## Stdev 0.0235
## Skewness -0.1070
## Kurtosis 5.9038
Plot showing portfolio returns and the dates of mass shootings (vertical dash lines) on function of time
gun_portfolio_returns_plot <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra))
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_line(aes(color="red"), size=1)
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + labs(title = "Returns of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + ylab("Ra") + ylim(-0.2, 0.2)
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_hline(yintercept = 0, color="black")
gun_portfolio_returns_plot
ggsave("89plot09B.png")
## Saving 7 x 5 in image
Bar plot showing the daily return of gun industry portfolio
gun_portfolio_returns_bars <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra))
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_bar(stat = "identity", fill=palette_light()[[1]])
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_smooth(method = "lm")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + labs(title = "Returns of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + ylab("Ra") + ylim(-0.05,0.05)
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_hline(yintercept = 0, color="black")
gun_portfolio_returns_bars
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 140 rows containing non-finite values (stat_smooth).
## Warning: Removed 140 rows containing missing values (position_stack).
ggsave("89plot10B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 140 rows containing non-finite values (stat_smooth).
## Removed 140 rows containing missing values (position_stack).
Plot showing portfolio drawdowns and mass shootings (vertical das lines) on function of time
gun_portfolio_DD_plot <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=gunDD))
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + geom_line(aes(color="red"), size=1)
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + labs(title = "DrawDowns of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + ylab("Draw Downs") + ylim(-0.5, 0)
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_DD_plot
ggsave("89plot11B.png")
## Saving 7 x 5 in image
Plot showing 5 day moving average and mass shooting events on function of time
gun_portfolio_RaAvg <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra_5day))
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_line(aes(color="red"), size=1)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + labs(title = "5 day moving average of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + ylab("Draw Downs")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_RaAvg
ggsave("89plot12B.png")
## Saving 7 x 5 in image
Plot showing 10 day portfolio momentun and mass shootign events on function of time
gun_portfolio_RaAvg <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Mom_10day))
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_line(aes(color="red"), size=1)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + labs(title = "10 day return momentum of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + ylab("Draw Downs")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_RaAvg
## Warning: Removed 10 row(s) containing missing values (geom_path).
ggsave("89plot13B.png")
## Saving 7 x 5 in image
## Warning: Removed 10 row(s) containing missing values (geom_path).
Cumulative return using tq_portfolio function
gun_co_portfolio_growth <- gun_co_returns %>% tq_portfolio(assets_col = symbol, returns_col = Ra, weights = gun_wts, col_rename = "investment.growth", wealth.index=TRUE)
gun_co_portfolio_growth$investment.growth <- gun_co_portfolio_growth$investment.growth-1
head(gun_co_portfolio_growth)
## # A tibble: 6 × 2
## date investment.growth
## <date> <dbl>
## 1 2009-10-01 0
## 2 2009-10-02 -0.0272
## 3 2009-10-05 -0.0350
## 4 2009-10-06 -0.0163
## 5 2009-10-07 -0.00777
## 6 2009-10-08 -0.00000128
Plot showing cumulative portfolio return and mass shooitng events on function of time
gun_portfolio_cum_return <- ggplot(data = gun_co_portfolio_growth, aes(x=date, y=investment.growth))
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_line(size=1, color=palette_light()[[1]])
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_smooth(method = "loess")
gun_portfolio_cum_return <- gun_portfolio_cum_return + theme(legend.position = "none")
## Vertical lines to show mass shootings
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)
gun_portfolio_cum_return <- gun_portfolio_cum_return + labs(title = "Compounded Return of US Gun Manufacture Portfolio",
subtitle = "Mass Shootings with seven or more total victims",
caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_cum_return <- gun_portfolio_cum_return + ylab("Ra")
gun_portfolio_cum_return <- gun_portfolio_cum_return + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_cum_return
## `geom_smooth()` using formula 'y ~ x'
ggsave("89plot14B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
Merging mass shooting data set and portfolio data set. Those observations where 20 day moving average is NA were eliminated. An indicator (dummy) variable for shooting was created
temp_mass_dates <- us_mass02 %>% filter(total_victims >= 7) %>% select("date", "fatalities", "injured", "total_victims")
merged_df <- merge(gun_co_portfolio_returns, temp_mass_dates, by.x = "date", by.y = "date", all.x = TRUE)
merged_df <- merged_df %>% replace_na(list(fatalities=0, injured=0, total_victims=0)) %>%
mutate(shooting = ifelse(total_victims>0,1,0))
merged_df2 <- merged_df %>% drop_na(Mom_20day) #%>% drop_na(Ra_10day)
tail(merged_df2, 10)
## date Ra SMA.2day SMA.15day gunDD gunDD_2day
## 3132 2022-06-09 4.580987e-03 -0.0041478253 0.004347252 -0.3335843 -0.3554429
## 3133 2022-06-10 -2.838113e-03 0.0008714373 0.004435712 -0.3354756 -0.3554188
## 3134 2022-06-13 -3.004751e-02 -0.0164428100 0.001333129 -0.3554429 -0.3522838
## 3135 2022-06-14 3.746473e-05 -0.0150050213 0.001744966 -0.3554188 -0.3804124
## 3136 2022-06-15 4.863690e-03 0.0024505774 0.001966289 -0.3522838 -0.3838896
## 3137 2022-06-16 -4.342748e-02 -0.0192818958 -0.004065449 -0.3804124 -0.3841058
## 3138 2022-06-17 -5.611991e-03 -0.0245197364 -0.005448874 -0.3838896 -0.3927233
## 3139 2022-06-21 -3.510373e-04 -0.0029815143 -0.005471049 -0.3841058 -0.3599289
## 3140 2022-06-22 -1.399185e-02 -0.0071714430 -0.007230260 -0.3927233 -0.3262850
## 3141 2022-06-23 5.400243e-02 0.0200052899 -0.003036847 -0.3599289 -0.3580488
## gunDD_3day gunDD_5day gunDD_10day gunDD_15day gunDD_20day Ra_2day
## 3132 -0.3554188 -0.3804124 -0.3262850 -0.3679977 -0.3731047 -0.0041478253
## 3133 -0.3522838 -0.3838896 -0.3580488 -0.3491225 -0.3785227 0.0008714373
## 3134 -0.3804124 -0.3841058 -0.3635919 -0.3760208 -0.3752875 -0.0164428100
## 3135 -0.3838896 -0.3927233 -0.3728398 -0.3685268 -0.3863849 -0.0150050213
## 3136 -0.3841058 -0.3599289 -0.3768972 -0.3757678 -0.3820448 0.0024505774
## 3137 -0.3927233 -0.3262850 -0.3679977 -0.3731047 -0.3950705 -0.0192818958
## 3138 -0.3599289 -0.3580488 -0.3491225 -0.3785227 -0.3872490 -0.0245197364
## 3139 -0.3262850 -0.3635919 -0.3760208 -0.3752875 -0.3773754 -0.0029815143
## 3140 -0.3580488 -0.3728398 -0.3685268 -0.3863849 -0.3727382 -0.0071714430
## 3141 -0.3635919 -0.3768972 -0.3757678 -0.3820448 -0.3760701 0.0200052899
## Ra_3day Ra_5day Ra_10day Mom_2day Mom_3day
## 3132 -4.185704e-05 -0.008620005 0.001043022 -0.034628495 -0.004543523
## 3133 -3.711254e-03 -0.003272175 -0.000754728 0.002875577 0.007701803
## 3134 -9.434878e-03 -0.006602238 -0.003757638 0.034911197 -0.013379974
## 3135 -1.094939e-02 -0.008228761 -0.004993524 -0.043464946 -0.005649456
## 3136 -8.382117e-03 -0.004680696 -0.003617277 -0.010475681 -0.005214727
## 3137 -1.284211e-02 -0.014282389 -0.011451197 0.043076444 0.029435633
## 3138 -1.472526e-02 -0.014837165 -0.009054670 -0.008379857 0.059614420
## 3139 -1.646350e-02 -0.008897871 -0.007750055 0.054353466 0.052913839
## 3140 -6.651626e-03 -0.011703734 -0.009966247 0.066554650 -0.033155442
## 3141 1.321985e-02 -0.001875986 -0.003278341 -0.101149720 -0.062637205
## Mom_5day Mom_10day Mom_15day Mom_20day fatalities injured
## 3132 -0.0480084690 0.04798181 0.009701543 -0.0003147059 0 0
## 3133 -0.0027738785 -0.04430918 0.032703896 -0.0058045291 0 0
## 3134 0.0296964700 0.02141273 -0.011278637 0.0352531855 0 0
## 3135 -0.0140293134 -0.01456875 0.011972425 -0.0178014742 0 0
## 3136 0.0491387385 -0.01133326 -0.016330538 0.0022093910 0 0
## 3137 0.0959902830 0.05771001 0.047693763 0.0223487799 0 0
## 3138 -0.0415352998 0.03547777 -0.003030651 0.0185414733 0 0
## 3139 -0.0082837389 -0.04097511 0.005556716 0.0164646274 0 0
## 3140 -0.0005394359 0.02600174 -0.003772161 0.0214397378 0 0
## 3141 -0.0604720002 -0.06546928 -0.046929348 -0.0593143249 0 0
## total_victims shooting
## 3132 0 0
## 3133 0 0
## 3134 0 0
## 3135 0 0
## 3136 0 0
## 3137 0 0
## 3138 0 0
## 3139 0 0
## 3140 0 0
## 3141 0 0
A correlation plot was set up in order to determine multicollinearity among variables
corrplot::corrplot(cor(merged_df2[,-1]), method="number", number.cex = 0.7) #addCoef.col=1,
ggsave("89plot15B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
The data set was split into train and test. Until December 31, 2019 was used as a training, after that was used for testing
merged_df_train <- merged_df2 %>% filter(date<="2019-12-31")
merged_df_test <- merged_df2 %>% filter(date > "2019-12-31")
tail(merged_df_test, 10)
## date Ra SMA.2day SMA.15day gunDD gunDD_2day
## 615 2022-06-09 4.580987e-03 -0.0041478253 0.004347252 -0.3335843 -0.3554429
## 616 2022-06-10 -2.838113e-03 0.0008714373 0.004435712 -0.3354756 -0.3554188
## 617 2022-06-13 -3.004751e-02 -0.0164428100 0.001333129 -0.3554429 -0.3522838
## 618 2022-06-14 3.746473e-05 -0.0150050213 0.001744966 -0.3554188 -0.3804124
## 619 2022-06-15 4.863690e-03 0.0024505774 0.001966289 -0.3522838 -0.3838896
## 620 2022-06-16 -4.342748e-02 -0.0192818958 -0.004065449 -0.3804124 -0.3841058
## 621 2022-06-17 -5.611991e-03 -0.0245197364 -0.005448874 -0.3838896 -0.3927233
## 622 2022-06-21 -3.510373e-04 -0.0029815143 -0.005471049 -0.3841058 -0.3599289
## 623 2022-06-22 -1.399185e-02 -0.0071714430 -0.007230260 -0.3927233 -0.3262850
## 624 2022-06-23 5.400243e-02 0.0200052899 -0.003036847 -0.3599289 -0.3580488
## gunDD_3day gunDD_5day gunDD_10day gunDD_15day gunDD_20day Ra_2day
## 615 -0.3554188 -0.3804124 -0.3262850 -0.3679977 -0.3731047 -0.0041478253
## 616 -0.3522838 -0.3838896 -0.3580488 -0.3491225 -0.3785227 0.0008714373
## 617 -0.3804124 -0.3841058 -0.3635919 -0.3760208 -0.3752875 -0.0164428100
## 618 -0.3838896 -0.3927233 -0.3728398 -0.3685268 -0.3863849 -0.0150050213
## 619 -0.3841058 -0.3599289 -0.3768972 -0.3757678 -0.3820448 0.0024505774
## 620 -0.3927233 -0.3262850 -0.3679977 -0.3731047 -0.3950705 -0.0192818958
## 621 -0.3599289 -0.3580488 -0.3491225 -0.3785227 -0.3872490 -0.0245197364
## 622 -0.3262850 -0.3635919 -0.3760208 -0.3752875 -0.3773754 -0.0029815143
## 623 -0.3580488 -0.3728398 -0.3685268 -0.3863849 -0.3727382 -0.0071714430
## 624 -0.3635919 -0.3768972 -0.3757678 -0.3820448 -0.3760701 0.0200052899
## Ra_3day Ra_5day Ra_10day Mom_2day Mom_3day
## 615 -4.185704e-05 -0.008620005 0.001043022 -0.034628495 -0.004543523
## 616 -3.711254e-03 -0.003272175 -0.000754728 0.002875577 0.007701803
## 617 -9.434878e-03 -0.006602238 -0.003757638 0.034911197 -0.013379974
## 618 -1.094939e-02 -0.008228761 -0.004993524 -0.043464946 -0.005649456
## 619 -8.382117e-03 -0.004680696 -0.003617277 -0.010475681 -0.005214727
## 620 -1.284211e-02 -0.014282389 -0.011451197 0.043076444 0.029435633
## 621 -1.472526e-02 -0.014837165 -0.009054670 -0.008379857 0.059614420
## 622 -1.646350e-02 -0.008897871 -0.007750055 0.054353466 0.052913839
## 623 -6.651626e-03 -0.011703734 -0.009966247 0.066554650 -0.033155442
## 624 1.321985e-02 -0.001875986 -0.003278341 -0.101149720 -0.062637205
## Mom_5day Mom_10day Mom_15day Mom_20day fatalities injured
## 615 -0.0480084690 0.04798181 0.009701543 -0.0003147059 0 0
## 616 -0.0027738785 -0.04430918 0.032703896 -0.0058045291 0 0
## 617 0.0296964700 0.02141273 -0.011278637 0.0352531855 0 0
## 618 -0.0140293134 -0.01456875 0.011972425 -0.0178014742 0 0
## 619 0.0491387385 -0.01133326 -0.016330538 0.0022093910 0 0
## 620 0.0959902830 0.05771001 0.047693763 0.0223487799 0 0
## 621 -0.0415352998 0.03547777 -0.003030651 0.0185414733 0 0
## 622 -0.0082837389 -0.04097511 0.005556716 0.0164646274 0 0
## 623 -0.0005394359 0.02600174 -0.003772161 0.0214397378 0 0
## 624 -0.0604720002 -0.06546928 -0.046929348 -0.0593143249 0 0
## total_victims shooting
## 615 0 0
## 616 0 0
## 617 0 0
## 618 0 0
## 619 0 0
## 620 0 0
## 621 0 0
## 622 0 0
## 623 0 0
## 624 0 0
The following lines are the series of models built for this analysis. Transformations were implemented as well
Model1 <- lm(data = merged_df_train, formula = SMA.2day ~ fatalities + injured)
summary(Model1)
##
## Call:
## lm(formula = SMA.2day ~ fatalities + injured, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.129373 -0.008098 0.000027 0.008098 0.086079
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.606e-04 3.226e-04 2.667 0.00769 **
## fatalities 3.220e-05 2.240e-04 0.144 0.88570
## injured 1.756e-05 4.001e-05 0.439 0.66083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01611 on 2514 degrees of freedom
## Multiple R-squared: 0.0002286, Adjusted R-squared: -0.0005668
## F-statistic: 0.2874 on 2 and 2514 DF, p-value: 0.7502
plot(Model1)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
Model5 <- lm(data = merged_df_train, formula = SMA.2day ~ total_victims)
summary(Model5)
##
## Call:
## lm(formula = SMA.2day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.129375 -0.008100 0.000025 0.008096 0.086077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.621e-04 3.214e-04 2.683 0.00735 **
## total_victims 1.934e-05 2.557e-05 0.756 0.44967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01611 on 2515 degrees of freedom
## Multiple R-squared: 0.0002273, Adjusted R-squared: -0.0001703
## F-statistic: 0.5717 on 1 and 2515 DF, p-value: 0.4497
plot(Model5)
Model2 <- lm(data = merged_df_train, formula = SMA.15day ~ fatalities + injured)
summary(Model2)
##
## Call:
## lm(formula = SMA.15day ~ fatalities + injured, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0219520 -0.0039910 0.0001739 0.0037688 0.0198369
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.870e-04 1.165e-04 7.614 3.73e-14 ***
## fatalities -1.879e-04 8.089e-05 -2.323 0.02027 *
## injured 4.036e-05 1.445e-05 2.793 0.00526 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005818 on 2514 degrees of freedom
## Multiple R-squared: 0.003211, Adjusted R-squared: 0.002418
## F-statistic: 4.049 on 2 and 2514 DF, p-value: 0.01756
par(mfrow=c(1,2))
plot(Model2)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
require(ggiraph)
## Loading required package: ggiraph
## Warning: package 'ggiraph' was built under R version 4.2.1
require(ggiraphExtra)
## Loading required package: ggiraphExtra
## Warning: package 'ggiraphExtra' was built under R version 4.2.1
require(plyr)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following object is masked from 'package:purrr':
##
## compact
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
#install.packages("devtools")
#devtools::install_github("cardiomoon/ggiraphExtra")
#install.packages("predict3d")
ggPredict(Model2,se=TRUE,interactive=TRUE)
Model6 <- lm(data = merged_df_train, formula = SMA.15day ~ total_victims)
summary(Model6)
##
## Call:
## lm(formula = SMA.15day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0219273 -0.0039855 0.0001545 0.0037848 0.0198616
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.624e-04 1.162e-04 7.421 1.58e-13 ***
## total_victims 1.264e-05 9.247e-06 1.366 0.172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005824 on 2515 degrees of freedom
## Multiple R-squared: 0.0007419, Adjusted R-squared: 0.0003446
## F-statistic: 1.867 on 1 and 2515 DF, p-value: 0.1719
plot(Model6)
DD1 <- lm(data = merged_df_train, formula = gunDD_2day ~ total_victims)
summary(DD1)
##
## Call:
## lm(formula = gunDD_2day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33012 -0.11365 0.01923 0.10485 0.22278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2175150 0.0027745 -78.397 <2e-16 ***
## total_victims -0.0002633 0.0002208 -1.193 0.233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1391 on 2515 degrees of freedom
## Multiple R-squared: 0.0005654, Adjusted R-squared: 0.000168
## F-statistic: 1.423 on 1 and 2515 DF, p-value: 0.2331
DD2 <- lm(data = merged_df_train, formula = gunDD_3day ~ total_victims)
summary(DD2)
##
## Call:
## lm(formula = gunDD_3day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33002 -0.11396 0.01924 0.10494 0.22180
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2176200 0.0027756 -78.404 <2e-16 ***
## total_victims -0.0002459 0.0002209 -1.113 0.266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1391 on 2515 degrees of freedom
## Multiple R-squared: 0.0004924, Adjusted R-squared: 9.498e-05
## F-statistic: 1.239 on 1 and 2515 DF, p-value: 0.2658
DD3 <- lm(data = merged_df_train, formula = gunDD_5day ~ total_victims)
summary(DD3)
##
## Call:
## lm(formula = gunDD_5day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32985 -0.11403 0.01925 0.10523 0.22367
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2177914 0.0027778 -78.405 <2e-16 ***
## total_victims -0.0002937 0.0002210 -1.329 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1392 on 2515 degrees of freedom
## Multiple R-squared: 0.0007016, Adjusted R-squared: 0.0003043
## F-statistic: 1.766 on 1 and 2515 DF, p-value: 0.184
DD4 <- lm(data = merged_df_train, formula = gunDD_10day ~ total_victims)
summary(DD4)
##
## Call:
## lm(formula = gunDD_10day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32939 -0.11425 0.01935 0.10557 0.21825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2182529 0.0027830 -78.425 <2e-16 ***
## total_victims -0.0003128 0.0002215 -1.413 0.158
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1395 on 2515 degrees of freedom
## Multiple R-squared: 0.0007928, Adjusted R-squared: 0.0003955
## F-statistic: 1.995 on 1 and 2515 DF, p-value: 0.1579
DD5 <- lm(data = merged_df_train, formula = gunDD_15day ~ total_victims)
summary(DD5)
##
## Call:
## lm(formula = gunDD_15day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32888 -0.11420 0.01863 0.10624 0.22320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2187592 0.0027882 -78.459 <2e-16 ***
## total_victims -0.0002610 0.0002219 -1.176 0.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1397 on 2515 degrees of freedom
## Multiple R-squared: 0.0005498, Adjusted R-squared: 0.0001524
## F-statistic: 1.383 on 1 and 2515 DF, p-value: 0.2396
DD6 <- lm(data = merged_df_train, formula = gunDD_20day ~ total_victims)
summary(DD6)
##
## Call:
## lm(formula = gunDD_20day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32846 -0.11512 0.01825 0.10651 0.21918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2191758 0.0027932 -78.467 <2e-16 ***
## total_victims -0.0002939 0.0002223 -1.322 0.186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.14 on 2515 degrees of freedom
## Multiple R-squared: 0.0006949, Adjusted R-squared: 0.0002976
## F-statistic: 1.749 on 1 and 2515 DF, p-value: 0.1861
mom1 <- lm(data = merged_df_train, formula = Mom_2day ~ total_victims)
summary(mom1)
##
## Call:
## lm(formula = Mom_2day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19071 -0.01610 -0.00056 0.01634 0.19255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.598e-05 6.384e-04 0.056 0.9551
## total_victims -1.147e-04 5.080e-05 -2.257 0.0241 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03199 on 2515 degrees of freedom
## Multiple R-squared: 0.002022, Adjusted R-squared: 0.001625
## F-statistic: 5.095 on 1 and 2515 DF, p-value: 0.02409
mom2 <- lm(data = merged_df_train, formula = Mom_3day ~ total_victims)
summary(mom2)
##
## Call:
## lm(formula = Mom_3day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.170502 -0.017780 -0.000561 0.017318 0.229605
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.330e-05 6.551e-04 0.020 0.984
## total_victims -4.235e-05 5.213e-05 -0.812 0.417
##
## Residual standard error: 0.03283 on 2515 degrees of freedom
## Multiple R-squared: 0.0002624, Adjusted R-squared: -0.0001351
## F-statistic: 0.66 on 1 and 2515 DF, p-value: 0.4166
Model3 <- lm(data = merged_df_train, formula = Mom_5day ~ fatalities + injured)
summary(Model3)
##
## Call:
## lm(formula = Mom_5day ~ fatalities + injured, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.169070 -0.016474 -0.000299 0.017039 0.191508
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.270e-05 6.530e-04 0.081 0.936
## fatalities -4.435e-05 4.534e-04 -0.098 0.922
## injured -1.295e-04 8.099e-05 -1.599 0.110
##
## Residual standard error: 0.03261 on 2514 degrees of freedom
## Multiple R-squared: 0.002112, Adjusted R-squared: 0.001318
## F-statistic: 2.661 on 2 and 2514 DF, p-value: 0.0701
plot(Model3)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
Model7 <- lm(data = merged_df_train, formula = Mom_5day ~ total_victims)
summary(Model7)
##
## Call:
## lm(formula = Mom_5day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.169079 -0.016358 -0.000308 0.017030 0.191499
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.189e-05 6.505e-04 0.095 0.9242
## total_victims -1.191e-04 5.177e-05 -2.301 0.0215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0326 on 2515 degrees of freedom
## Multiple R-squared: 0.002101, Adjusted R-squared: 0.001704
## F-statistic: 5.296 on 1 and 2515 DF, p-value: 0.02146
par(mfrow=c(2,2))
plot(Model7)
ggPredict(Model7,se=TRUE,interactive=TRUE)
mom4 <- lm(data = merged_df_train, formula = Mom_10day ~ total_victims)
summary(mom4)
##
## Call:
## lm(formula = Mom_10day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.174663 -0.017609 0.000211 0.017679 0.163871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.094e-05 6.423e-04 0.033 0.974
## total_victims -4.767e-05 5.111e-05 -0.933 0.351
##
## Residual standard error: 0.03219 on 2515 degrees of freedom
## Multiple R-squared: 0.0003457, Adjusted R-squared: -5.176e-05
## F-statistic: 0.8698 on 1 and 2515 DF, p-value: 0.3511
mom5 <- lm(data = merged_df_train, formula = Mom_15day ~ total_victims )
summary(mom5)
##
## Call:
## lm(formula = Mom_15day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.171799 -0.017124 0.000485 0.017619 0.191330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.578e-05 6.491e-04 0.071 0.9438
## total_victims -9.134e-05 5.165e-05 -1.768 0.0771 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03253 on 2515 degrees of freedom
## Multiple R-squared: 0.001242, Adjusted R-squared: 0.0008448
## F-statistic: 3.127 on 1 and 2515 DF, p-value: 0.07711
Model4 <- lm(data = merged_df_train, formula = Mom_20day ~ fatalities + injured)
summary(Model4)
##
## Call:
## lm(formula = Mom_20day ~ fatalities + injured, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.181020 -0.017133 -0.000324 0.016720 0.206396
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.075e-04 6.569e-04 0.164 0.870
## fatalities -4.548e-04 4.561e-04 -0.997 0.319
## injured -6.589e-05 8.147e-05 -0.809 0.419
##
## Residual standard error: 0.0328 on 2514 degrees of freedom
## Multiple R-squared: 0.002098, Adjusted R-squared: 0.001304
## F-statistic: 2.643 on 2 and 2514 DF, p-value: 0.07134
plot(Model4)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
Model8 <- lm(data = merged_df_train, formula = Mom_20day ~ total_victims)
summary(Model8)
##
## Call:
## lm(formula = Mom_20day ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.180978 -0.017098 -0.000289 0.016761 0.206438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.552e-05 6.545e-04 0.100 0.9203
## total_victims -1.131e-04 5.208e-05 -2.172 0.0299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0328 on 2515 degrees of freedom
## Multiple R-squared: 0.001873, Adjusted R-squared: 0.001476
## F-statistic: 4.718 on 1 and 2515 DF, p-value: 0.02993
plot(Model8)
Model8a <- lm(data = merged_df_train, formula = Mom_20day ~ log(total_victims+1))
summary(Model8a)
##
## Call:
## lm(formula = Mom_20day ~ log(total_victims + 1), data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.181023 -0.017137 -0.000328 0.016766 0.206392
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0001111 0.0006592 0.169 0.866
## log(total_victims + 1) -0.0022828 0.0016907 -1.350 0.177
##
## Residual standard error: 0.03282 on 2515 degrees of freedom
## Multiple R-squared: 0.0007243, Adjusted R-squared: 0.000327
## F-statistic: 1.823 on 1 and 2515 DF, p-value: 0.1771
c <- abs(min(merged_df_train$Mom_20day))
c
## [1] 0.180912
merged_df_train$log_mom_20day_c <- log(merged_df_train$Mom_20day+c+1)
Model8b <- lm(data = merged_df_train, formula = log_mom_20day_c ~ total_victims)
summary(Model8b)
##
## Call:
## lm(formula = log_mom_20day_c ~ total_victims, data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.165957 -0.014198 0.000141 0.014478 0.161486
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.660e-01 5.556e-04 298.716 <2e-16 ***
## total_victims -9.839e-05 4.421e-05 -2.225 0.0261 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02784 on 2515 degrees of freedom
## Multiple R-squared: 0.001965, Adjusted R-squared: 0.001569
## F-statistic: 4.953 on 1 and 2515 DF, p-value: 0.02614
par(mfrow=c(1,2))
plot(Model8b)
ggPredict(Model8b, se=TRUE,interactive=TRUE)
c <- abs(min(merged_df_train$Mom_20day))
c
## [1] 0.180912
Model8c <- lm(data = merged_df_train, formula = log(Mom_20day+c+1) ~ log(total_victims+1))
summary(Model8c)
##
## Call:
## lm(formula = log(Mom_20day + c + 1) ~ log(total_victims + 1),
## data = merged_df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.166000 -0.014235 0.000104 0.014478 0.161443
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1659997 0.0005596 296.62 <2e-16 ***
## log(total_victims + 1) -0.0020517 0.0014352 -1.43 0.153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02786 on 2515 degrees of freedom
## Multiple R-squared: 0.0008119, Adjusted R-squared: 0.0004146
## F-statistic: 2.044 on 1 and 2515 DF, p-value: 0.153
library(forecast)
predict_01 <- predict(Model2, newdata = merged_df_test)
head(predict_01,10)
## 1 2 3 4 5 6
## 0.0008870347 0.0008870347 0.0008870347 0.0008870347 0.0008870347 0.0008870347
## 7 8 9 10
## 0.0008870347 0.0008870347 0.0008870347 0.0008870347
accuracy_01 <- accuracy(predict_01, merged_df_test$SMA.15day)
accuracy_01
## ME RMSE MAE MPE MAPE
## Test set 0.000321654 0.005971153 0.004634279 82.64277 140.5359
predict_02 <- predict(Model7, newdata = merged_df_test)
head(predict_02,10)
## 1 2 3 4 5 6
## 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05
## 7 8 9 10
## 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05
accuracy_02 <- accuracy(predict_02, merged_df_test$Mom_5day)
accuracy_02
## ME RMSE MAE MPE MAPE
## Test set -0.0001480064 0.03681865 0.02652235 93.54564 106.0142
gun_co_portfolio_returns1 <- gun_co_portfolio_returns
gun_co_portfolio_returns1 <- xts(gun_co_portfolio_returns1[,-1], order.by = gun_co_portfolio_returns$date,)
#head(gun_co_portfolio_returns1)
gun_comp_return <- Return.cumulative(gun_co_portfolio_returns1$Ra, geometric = TRUE)
gun_comp_return
## Ra
## Cumulative Return 6.729115
chart.CumReturns(gun_co_portfolio_returns1$Ra, wealth.index = FALSE, geometric = TRUE)
library(data.table)
gun_DD_plot <- chart.Drawdown(gun_co_portfolio_returns1$Ra)
gun_DD_plot
gun_DDtable <- table.Drawdowns(gun_co_portfolio_returns1$Ra, top = 7, digits = 4)
gun_DDtable
## From Trough To Depth Length To Trough Recovery
## 1 2014-01-17 2014-12-16 2016-02-26 -0.5476 531 231 300
## 2 2016-03-23 2019-09-03 2020-07-02 -0.5389 1078 868 210
## 3 2021-07-01 2022-07-18 <NA> -0.3951 268 263 NA
## 4 2012-05-03 2012-06-07 2012-11-23 -0.3630 141 25 116
## 5 2012-11-30 2012-12-18 2013-09-11 -0.3122 196 13 183
## 6 2011-08-31 2011-10-03 2011-12-21 -0.3036 79 23 56
## 7 2020-08-11 2020-11-24 2021-06-01 -0.2743 203 75 128